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ABSTRACT 

We report on NuSTAR , XMM-Newton and Swift observations of the gamma-ray binary 
1FGL J1018.6—5856. We measure the orbital period to be 16.544 ± 0.008 days using Swift data 
spanning 1900 days. The orbital per iod is different fr om the 2011 gamma-ray measurement which was 
used in the previous X-ray study of I An et ahl (I2013T) using ~400days of Swift data, but is consistent 
with a new gamma-ray solution reported in 2014. The light curve folded on the new period is qualita¬ 
tively similar to that reported previously, having a spike at phase 0 and broad sinusoidal modulation. 

The X-ray flux enhancement at phase 0 occurs more regularly in time than was previously suggested. 

A spiky structure at this phase seems to be a persistent feature, although there is some variability. 
Furthermore, we find that the source flux clearly correlates with the spectral hardness throughout all 
orbital phases, and that the broadband X-ray spectra measured with NuSTAR , XMM-Newton , and 
Swift are well fit with an unbroken power-law model. This spectrum suggests that the system may 
not be accretion-powered. 

Subject headings: binaries: close — gamma rays: stars — X-rays: binaries — stars: individual 
(1FGL J1018.6—5856) 


1. INTRODUCTION 

Gamma-ray binaries are systems composed of a mas¬ 
sive star and a compact object and from which persistent 
GeV and/or TeV gamma-ray emission is detected and 
dominates the overall non-thermal spectrum. They emit 
across the electrom agnetic spectru m from the radio to 
TeV gamma ray (see lMirabel II2012I . for a review) . There 
are o nly five gamma-ray binaries known to date dDubus I 
EqU, and only for one source has the compact obje ct 
been identified (PSR B1259—63: lJohnston et al.lll992H . 

Since most of the energy output of a gamma-ray bi¬ 
nary is in the gamma-ray band, current theoretical stud¬ 
ies focus on explaining the high energy emission proper¬ 
ties. The gamma-ray emission models can be categorized 
into two classes: microquasar models le.g.. lRomero et all 
l2003HBosch-Ramon &; Paredesll2004lf and pulsar models 

(e.g.,rfhvani et al. III994HSierpowska-Bartoski fe Torres I 

I2008T) . In the microquasar model, relativistic elec¬ 
trons in a jet generated close to the compact ob¬ 
ject Compton-upscatter the synchrotron emission of 
the jet itself and/or the stellar UV p hotons (e.g., 
iKaufman-Bernado et a.1.1 l2002fc iBosch-Ramon &; Pareded 


H®, or relativistic hadrons collide with background 
nucle i creating pions that decay (e.g., iRomero et al.1 
1200311 . producing gamma rays. In the pulsar model, 
pulsar wind particles are accelerated in the pulsar 
wind/stellar wind shock, and Compton-upscatter stel¬ 
lar photons to produce the observed gamma rays (e.g. , 
Tavani et al. 111994 [Tavani fc Arons1ll997tlDubus 112006! 

Sieroowska-Bartoski fc Torres 1120081b 
Non-thermal X-ray emission in gamma-ray binaries 
is thought to be produced by the electrons which 
are accelerated in the pulsar wind/stellar wind shock 
(e.g. iTavani fe Arons I Il997t IDubus 1 120061) or in rela¬ 
tivistic jets formed close to the c ompact object (e.g., 
iBosch-Ramon fc Khangulvanl l2009 f. The models pre¬ 
dict varying X-ray fluxes and spectra depending on 
the properties of the shock, which are determined 
by the thrust of the winds and the orbita l geom¬ 
etry of the binary system (e.g., iKaspi et all 119951) . 
or on the jet dynamics and cooling timescale (e.g . , 
IDubus et all l201(t IBosch-Ramon fc Khangulvanl [2009 1. 
Hence, X-ray measurements can be used for constrain¬ 
ing the orbital parameters and understanding the na¬ 
ture of the physical processes in gamma-ray binaries (see 
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also iChernvakova et al. I 120061 : iTakahashi et al.1 120091 : 
iTakata et al. 1120 HF 

The gamma-ray binary 1FG L J1018.6—5856 was dis - 
covered with Fermi in 2011. lAckermann et all 620121 1 
found modulation in the radio to gamma-ray bands with 
a period of 16.58 ± 0.02 days, identifying the source as a 
gamma-ray binary. They further identified the compan¬ 
ion star to be an 06V((f)) star. Soon after the d iscovery 
subsequent broadband studies were carried out (|Li et al.l 
1 201 l as lAbramowski ct al~ll2012fc An ct al.ll2013D in order 
to better characterize the source properties, but in no 
case were they able to identify the nature of the compact 
object. 

X-ray properties of the gamma-ray binary 
1FGL J1018.6 —5856 were measured in detail with 
Swift. lAn ct ah (120131) showed that the X-ray flux 
peak seen at phase 0 (gamma-ray maximum) by 
lAckermann et all (120121 ) seems not to be a persistent 
feature and instead sho ws a relativ ely large orbit-to-orbit 
variation. Furthermore, I An et all (120131) found evidence 
of a correlation between flux and spectral hardness in 
the X-ray b and. _ 

Recently, iColev et al. : (|2014f ) refined the gamma-ray 
period using Fermi observations with a longer baseline, 
and found the period to be 16.531 ± 0.006 days. Since 
this is slightly different from the value (16.58± 0.02days) 
used f or the previous X-ray study carried out bv lAn et al.l 
( 2013 1. the X-ray results need to be refined using the new 
gamma-ray period. The baseline of the X-ray observa¬ 
tions is long (5 years), and thus phases of later observa¬ 
tions may change significantly. 

Important questions to be addressed for gamma-ray 
binaries are: what is the nature of the compact objec t 
(known only for PSR B1259—63. 1.Tohnston et aI1ll992f ). 
and what is the physical emission mechanism. If the 
source is powered by accretion, a complex continuum 
spectrum is expected whether the compact object is a 
neutron star or a black hole. Hence, accurate mea¬ 
surement of the spectrum will help us identify the com¬ 
pact object. Furthermore, searching for a spectral turn¬ 
over in the hard X -ray band (e.g., iGrove et al.l 119991 : 
iGoburn et al. I [20021) and/or spectral lines often seen in 
high-mass X-ray binaries (HMXBs) may also provide 
clues about the emission mechanism of the source. 

In this paper, we measure X-ray properties of the 
gamma-ray binary 1FGL J1018.6—5856 more accurately 
than before using new observations taken with NuSTAR, 
Swift and with archival XMM-Newton observations. In 
Section [2j we describe the observations we used in this 
paper. We show data analysis and the results in Sec¬ 
tion [3j We then discuss our findings in Section [4j and 
conclude in Section [5] 

2. OBSERVATIONS 

We observed the gamma- ray binary 
1FGL J1018.6—5856 with NuSTAR (lHarrison et al.1 
I2013T) four times between 2014 June 4 and December 
1 with exposures of ~20ks for each observation. The 
total exposure was 90 ks. Soft X-ray band below 
and overlapping with the NuSTAR band (3-79 keV) 
was covered with Swift observations and two archival 
XMM-Newton observations (see Table [£]). The total 
exposure of the 71 Swift observations was 169 ks, and 
each exposure was relatively short. 


The NuSTAR observations were processed with the 
standard pipeline tools nupipeline and nuproducts of 
nustardas 1.4.1 integrated in HEASOFT 6.16. We used 
NuSTAR CALDB version 20140414 and applied the stan¬ 
dard filters!]] In order to process the Swift data, we used 
the xrtpipeline tool along with HEASA RC remote 
C ALDE0 and standard filters dCapalbi et al.l2005l) . Note 
that the source was not clearly detected in some Swift 
observations, and that the Swift obser vations taken until 
MJD 55984 were rep orted previously (|Ackermann et all 
l2012t lAn et al~ll20 1 .'ill . The XMM-Newton data were pro¬ 
cessed with epproc and emproc in Science Analysis Sys¬ 
tem (SAS) 14.0. (0 using standard filters. 

3. DATA ANALYSIS AND RESULTS 
3.1. Timing Analysis 

Detection of pulsations in gamma-ray binaries can be 
difficult for several reasons, such as the possibilities of an 
unfavorable emission geometry, absorption of soft X-rays 
by the wind, or a large background due to non-thermal 
unpulsed emission. Even in a favorable situation where 
the above effects are minimal, the Doppler effect due to 
binary motion can blur the pul se s ig nal i f the orbit is 
tight. For 1FGL J1018.6-5856. lAn et al.l (120131) showed 
that the Doppler broadening is not a concern for a 20-ks 
observation assuming a circular orbit with an inclina¬ 
tion of 30°. We therefore attempt to search for the pul¬ 
sation. Event arrival times measured at the spacecraft 
were transformed into those at the solar system barycen- 
ter with barycorr for the NuSTAR and barycen for the 
XMM-Newton data. We did not search the Swift data 
because of the paucity of counts in individual Swift ob¬ 
servations. 

For the NuSTAR data, we produced an event list for 
each observation in the 3-20 keV band using a circular 
aperture with R = 30". We performed the timing anal¬ 
ysis with the data from each NuSTAR focal plane mod- 
u 1(0 as well as with the combined dataset. Above 20 keV 
background dominates, and hence we adopt that as the 
high end of our band. Note that the results below do not 
depend strongly on the exact energy range or the aper¬ 
ture size. We folded the event time series to test periods 
betw een 10~ 4 — 10 3 s, and calculated Z\ (jBuccheri et al.1 
I1983D . We find that Z\ is fairly large for some test pe¬ 
riods. However, we find that the large Z\ seen in one 
observation is not reproduced in the others. We further 
verified that the large Z\ values are not significant. Note 
that the measured Z\ distribution does not follow a x 2 
distribution, but has a long tail, and thus we used a func¬ 
tional distribution obtained by fitting the measured Z/ 
distribution in order to estimate the significance. We per¬ 
formed the same study for the XMM-Newton/ PN data 
in the 0.5-2 keV and 0.5-10 keV bands, and did not find 
any significant pulsations. Assuming the pulse profile is 
a sine function with a period in the range of 0.1-1 s, we 
estimate the 90% upper limit for the pulse fractions to 
be 47% and 6% in the 3-20 keV and 0.5-10 keV bands, 


1 See http://heasarc.gsfc.nasa.gov/docs/nustar/analysis/mistar_swguide.pdf 
for more details 

2 http://heasarc.nasa.gov/docs/heasarc/caldb/caldb_remote_ac 
cess.html 

3 http://xmm.esac.esa.int/sas/ 

4 NuSTAR has two focal plane modules, FPMA and FPMB. 
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Table 1 

Summary of observations used in this work and results of spectral analysis 


Observatory 

Obs. ID 

Date 

(MJD) 

Exposure 

(ks) 

0 

JVh 

(10 22 cm —2 ) 

r 

F 3 -10 keV a 

( erg s —1 cm -2 ) 

Mode 

XMM 

0604700101 

55066 

20/12 

0.6 

0.65(5) 

1.65(7) 

5.1 ± 0.2 x 10 -13 

FW/FW b 

XMM 

0694390101 

56302 

104/73 

0.3 

0.72(2) 

1.57(2) 

1.09 ± 0.01 X 10 -12 

FW/SW b 

Swift 

00031912001- 

00090191001 

55103- 

56992 

169 

o 

1 

o 

d 

0.72 c 

1.2-1.8 

0.34 - 1.2 x 10~ 12 

PC 

NuSTAR 

30002020002 

56812 

22 

0.2 

0.72 c 

1.67(10) 

5.7 ± 0.5 X 10 -13 


NuSTAR 

30002020004 

56862 

23 

0.2 

0.72 c 

1.77(10) 

6.1 ± 0.5 X 10~ 13 


NuSTAR 

30002020006 

56911 

25 

0.2 

0.72° 

1.64(8) 

7.8 ± 0.5 X 10 -13 


NuSTAR 

30002020008 

56992 

21 

0.0 

0.72 c 

1.41(7) 

1.11 ± 0.06 X 10~ 12 



a Absorption-corrected flux. 

b For MOSl,2/PN. FW: Full window. SW: small window. 
c /V H was frozen for the Swift and NuSTAR data fit. 

respectively. 

Next, we refine the X-ray measurement of the orbital 
period by using a longer baseline using the Swift data 
over a longer time period than the previous work. Note 
that we did not use the data taken with XMM-Newton 
or NuSTAR because their count rate measurements can¬ 
not be_directlvcompared to those of Swift. As was done 
by I An et al.l (f2013h . we use epoch folding dLeahv 1119871) 
because of the unequal exposures of the observations. 
In the Swift observations, we extracted source and back¬ 
ground events in the 0.5-10 keV band within a 30" radius 
circle, and an annular region with inner radius 50" and 
outer radius 100", respectively. We then folded the event 
time series at test periods around P or b = 16.531 days 
dColev et ah II2014L producing a light curve with 16 bins. 
We used the sa me epoch for phase 0 as that used in th e 
previous studies (lAckcrmann et aO2012t lAn et ahll2013f ) . 
We calculated y 2 of the light curve for each trial period, 
and fo llowed the fitting technique as described in lLcahvl 
(Il987t ). Note, however, that we modeled the underly¬ 
ing continuum using a power-l aw fu n ction instead of the 
constant model employed by iLeahv I (I1987T ) , because the 
y 2 of the folded light curve is rising towards short pe¬ 
riods (see Figure [lj. The best-fit continuum model is 
Xcont ~ P 1 ': the exact value of the power-law index 
varies between 0.9 and 1.1 depending on the fit range. 

We find that the best-fit orbital period varies between 
16.538 and 16.55 days, depending on the number of bins, 
and the search step or search range. We varied the 
number of bins between 8 and 18 to ensure that the 
light curve is resolved and the y 2 statistic is applica¬ 
ble, and the search step between 0.001 and 0.005 days, 
smaller than the uncertainty in the P 0 rb measurement 
(AP or b < 0.01 days). The search range was varied be¬ 
tween ±1 day and ±15 days. We find that the variations 
are within ler of the measurements. The resulting period 
is Porb = 16.544 ± 0.008 days. We show the folded light 
curve in Figure [2ji. 

We also used the second harmonic to measure the or¬ 
bital period, since it can be measured with better pre¬ 
cision, and found that the measurement is more stable, 
varying only 6 x 10~ 4 days as a function of the number 
of bins, search step or the search range. The result is 
Porb = 16.543 ±0.004 days. Our result is consistent with 
the Fermi measure d value (P or b = 16.531 ± 0.006 days, 
iColcv et al. I I2014T) with a null hypothesis probability 
p = 0.09. 
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Figure 1. Chi-squared vs. frequency for the 0.5-10 keV Swift X- 
ray data. The search step is 0.01 day and, the number of phase 
bins is 16. Frequencies for the first three harmonics are denoted as 
blue dotted lines and the best-fit function is shown in red. 

3.2. Spectral Analysis 

Although Ah towards the source has been measured 
previously, the uncertainty was relatively large. Since 
there is one more XMM-Newton observation (Obs. ID 
069439 0101, Table UT ) taken after the previous X-ray 
study (lAn et al.ll2013f) . we can determine Ah more pre¬ 
cisely using the XMM-Newton observations. We ex¬ 
tracted the source spectrum from a circle with R = 16" 
(Obs. ID 0604700101) or R = 24" (Obs. ID 0694390101), 
and background spectra from a circle with R = 32" in 
a source-free region ~200" vertically upwards along the 
detector column from the source. Note that we used dif¬ 
ferent source extraction regions because of differences in 
exposure times. 

Since it has been suggested that the spectral hardness 
varies orbitally, we used different spectral slopes for ob¬ 
servations taken at different phases (Figure 0. Thus, we 
fit the two XMM-Newton spectra separately allowing all 
the fit parameters to vary. We grouped the spectra to 
have 20 counts per bin, and fit them with an absorbed 
power-law model with th e angr abundance in XSPEC 
(1 Anders & Clrevessell 198911 using y 2 statistics or l statis¬ 
tics dLoredo II1992H . The two methods yield consistent 
results. The best-fit Ah values for the observations are 
statistically consistent with each other (Table [0. Best- 
fit Ah values ob tained with a d ifferent abundanc e model 
(wilm in XSPEC; IWillms. Allen fc McCray 1120001) for the 
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Figure 2. X-ray light curve and the best-fit spectral parameters for a power-law model. Only Swift data are shown in panel (a) and 
measurements made with all three instruments are shown in panels (b) and (c). a) 0.5—lOkeV light curve as measured with Swift. Black 
dotted line shows the average count rate, black solid line is for the average source count ra te, and the blac k dashed line shows the average 
background co unt rate. Blue data points are measurements reported by (before 55985 MJD fXn et al.lf2013l) . and green data points are new 
measurements JColev et al. I [20141 . and this work). Cyan triangle denotes the “high-flux state”, two observations which had significantly 
larger count rates than the others at the same phase, and red diamonds are for time periods in which NuSTAR observations were made. 
b) 3-10 keV flux corrected for interstellar absorption. Data points for Swift , NuSTAR , and the XMM-Newton measurements are denoted 
in cross, triangle, and diamond, respectively, c) the best-fit photon index. Same symbols as in (b) are used in (c). 



Figure 3. Broadband X-ray spectra obtained with XMM-Newton , NuSTAR and Swift. Left: XMM-Newton spectra obtained at phases 
0.3 and 0.6. The harder and brighter spectra for phase 0.3 measured with PN, MOS1 and MOS2 are colored in blue, cyan and magenta, 
respectively and the spectra for phase 0.6 are colored in black, red and green (see Figures [2j) and c). Right: NuSTAR and Swift spectra 
for phase 0 without the high-flux state. Note that the Swift spectrum (green) is obtained by combining Swift observations taken at phase 
0 excluding the high-flux state. Black and red data points show the NuSTAR spectra measured with FPMA and FPMB, respectively. The 
best-fit power-law model is shown in solid lines in each panel. 


two spectra are still consistent with each other. There¬ 
fore, we use a common Nu value and find that a power- 
law model successfully explains the data (Figure |31 left). 
The best-fit value is Nu = 7.2 ± 0.2 x 10 21 cm -2 , and 
we use this value throughout this paper. Note that using 
jV h = 7.2 ± 0.2 x 10 21 cm -2 does not change the other 
spectral parameters in Table [T] significantly. We find that 
using the wilm abundance model changes the best-fit Nu 
values (to 0.93±0.08xl0 22 cm -2 , 1.03±0.02xl0 22 cm -2 , 
and 1.02 ± 0.02 x 10 22 cm' 2 for Obs. IDs 0604700101, 
0694390101 and combined, respectively), but the other 
spectral parameters do not change significantly. We note 
that the source count rates were less than 0.03-0.08 cps 
for MOS1/2, and 0.1-0.3 cps for PN, and hence pile-up 
is not a concern^ 

5 http: //xmm2.esac.esa.int/docs/documents / C AL-TN-0200-1- 


For the NuSTAR data, we extracted source and back¬ 
ground events from circular regions with R = 30" and 
R = 45", respectively. Backgrounds were extracted in 
the same detector chip as the source, offset ~4' from the 
source region. The source was detected above the back¬ 
ground up to 20-30 keV. We grouped the spectra to have 
a minimum of 20 counts per spectral bin, and used % 2 
statistics and l statistics; they provide consistent results. 

We jointly fit the data with a power-law model having 
different photon indexes for different orbital phases, and 
found that the best-fit parameters are T = 1.69±0.05 and 
F3-10 keV = 5.7 ± 0.4 x 10 -13 erg cm -2 s _1 for (j) = 0.2, 
and T = 1.41 ± 0.07 and ^3-10 keV = 1-11 ± 0.06 x 
10~ 12 erg cm -2 s _1 for (j> = 0. We find that a power- 
law model with a constant photon index across orbital 
phases can also explain the data if we let the flux vary 
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between phases. However, a model with separate spec¬ 
tral indices for different phases provides a significantly 
better fit than one with constant phase-independent in¬ 
dex throughout the phases, having an F-test probability 
that the improvement is just due to statistical chance of 
2 x 1CP 3 . Using separate power-law indexes for the three 
observations taken at cj> = 0.2 does not improve the fit. 
Furthermore, individual fits of the observations suggest 
that the photon index is statistically the same among the 
observations taken at (f> = 0.2 and that the photon index 
at phase 0 is different from that at phase 0.2. We show 
the NuSTAR spectra in Figure [3] and the fit results in 
Figures [2] and 01 

For the Swift data, we extract the spectra using the 
same regions as for the timing analysis (Section 13.11) . 
The center of the source extraction circle was deter¬ 
mined for each observation separately. Since t he source 
spect ral properties vary with orbital phase (lAn et alj 
mm . we performed phase-resolved spectroscopy. We 
folded the observations using the new timing solution 
we found in Section 13.11 and merged the data in each 
phase bin, for a total of twelve phase bins. We fur¬ 
ther produced two spectra for phase 0, one for the high- 
flux state and another for the rest of the observations 
taken at that phase, and eleven spectra for the other 
phases hence producing a total of thirteen spectra. We 
grouped the data to have 1 count per energy bin because 
of the paucity of counts in some phase bins, and used 
l statistics. For the phases that have enough counts, 
we also tried to fit the spectra using \ 2 statistics, af¬ 
ter grouping to have more than 20 counts per energy 
bin, and found that the results are consistent with those 
obtained using l statistics. We then fit all 13 spectra 
jointly with an absorbed power-law model with a com¬ 
mon Nu (frozen at the XMM-Newton-measmed value of 
7.2 ± 0.2 x 10 21 cm -2 ) throughout the observations but 
a separate photon index and flux for each spectrum. We 
find that the power-law model explains the data with 
photon indices of T = 1.2 — 1.8 and 3-10 keV fluxes of 
i< 3 _io kev = 0.34-2.9 x 10~ 12 erg cm -2 s -1 , where the 
maximum flux was for the high-flux state (see Figures [5] 
b and c). 

We also tried to determine Ah at each orbital phase 
using the Swift data. However, we were able to con¬ 
strain the fit parameters reasonably only for phase 0 
(without the high-flux state). We fit the spectrum for 
phase 0 (excluding the high-flux state) with a power- 
law model, and find that the best-fit parameters are 
A H = 7.7 ± 1.2 x 10 21 cm" 2 , T = 1.4 ± 0.1, and 
A 3 _io kev = 1.5 ± 0.1 x 10 -12 erg cm -2 s -1 . This Ah 
value at phase 0 is consistent with those obtained using 
the XMM-Newton data for phases 0.3 and 0.6 above, sug¬ 
gesting that IVh does not strongly vary as a function of 
orbital phase. Although we cannot clearly rule out or¬ 
bital variation of IVh, a —10% variation of Ah does not 
significantly change the Swift results. 

Accretion-powered neutron star HMXBs often show 
spectral features such as e mission lines or an ex ponential 
cutoff in the X-ray band dCoburn et al. I [20021) . We find 
that the X-ray spectrum of 1FGL J1018.6—5856 is well 
described with a power-law model without requiring any 
additional features (e.g., Figure [3]). For example, fitting 
the spectrum for phase 0 with a cutoff power-law model 


(pow*highecut in XSPEC) does not improve the fit, and 
the best-fit parameters are not constrained. We further 
changed the spectral grouping in order to have the spec¬ 
tra cover a broader energy range, and to see if a cutoff 
is required at higher energy. Specifically, we grouped the 
NuSTAR spectra to have more than 15 counts per energy 
bin, covering the 3-70 keV band. We fit the spectra with 
a power-law model and a cutoff power-law model, and 
found the same results as above; no cutoff is required in 
the fit. 

We performed additional analysis to determine the 
lower limit for the cutoff energy (A cu t 0 ff of the highecut 
model). However, it is not possible to set a meaning¬ 
ful lower limit for -E^uto® without constraining the e- 
folding energy ( Ef of the highecut model). We there¬ 
fore limit Ef bewteen 6 keV and 12 keV, values obtained 
for a sample of accret ion-powered neutron star HMXBs 
dCoburn et al. I 120021 1. and found that the 90% lower 
limit for A cuto ff is 39 keV and 34 keV for Ef of 6 keV and 
12keV, respectively. Note that some accretion-powered 
black hole binaries are known to have the cutoff energy 
above 70keV and that our data are not sensitive to such 
high energy cutoff. 

3.3. Spectral variability 

The spectral hardness varies with orbital phase (Fig¬ 
ure [2t) and flux (Figure [4}. Figure 0] shows an ap¬ 
parent correlation between flux and spectral hardness. 
We fit the apparent correlation with a constant func¬ 
tion and found that it does not provide an acceptable 
fit (x 2 /dof=55/18). We therefore added a linear slope to 
the constant function and find that the linear fit explains 
the data well (x 2 /dof=17/17). The measured slope is 
—0.28 ± 0.04 (pe r 10~ 12 erg cm~ 2 s -1 ), consistent with 
that reported bv lAn et~ (2013). The best-fit function 
is sh own in Figure [4] 

lAn et all (|2013f) suggested that there is evidence for 
a correlation between X-ray flux and spectral hardness. 
With the new and larger dataset, it is clear that the two 
quantities are correlated. For example, we find that the 
Spearman’s rank order correlation coefficient is —0.84 
and the significance is 5.6cr. We further verified that 
the X-ray flux and the photon index vary orbitally using 
the x 2 test, which resulted in p < 10~ 10 and p — 10 -5 , 
respectively. Note that whether or not we include the 
high-flux data point in the correlation calculation does 
not significantly change the result. 

Although the significance for the correlation is high, 
uncertainties in the measurements are significant (see 
Figure 0]) and need to be considered for the significance 
calculation. In order to do so, we performed simulations. 
Note that the photon index and the flux are correlated in 
the spectral fit, and one needs to take into account the 
covariance. We do this by using the covariance matri¬ 
ces in the simulation as was done by lAn et all (12013T) . 
In 100,000 simulations, a non-negative correlation oc¬ 
curred 316 times, which suggests that the significance 
of the negative correlation is —99.7%. We also carried 
out simulations for the linear correlation, and measured 
the confidence level of the negative linear correlation to 
be -99.9%. 

We also checked for short-term variability (—lOks) 
using the longer XMM-Newton observation (Obs. ID 
0694390101) because it has the best statistics. We cal- 
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Flux (10 12 erg s 1 cm 2 ) 

Figure 4. Photon index vs. 3-10 keV flux. The dashed blue line 
shows the best-fit linear function. 

culate the count rate and hardness ratio (ratio of count 
rates in two energy bands; e.g, C 3 -iokev/Co. 3 - 3 kev) on 
various time scales and energy bands. We find variabili¬ 
ties of ~20% and ~10% for the count rate and hardness 
ratio, respectively, but no correlation between them. 

4. DISCUSSION 

Our analysis of the NuSTAR , XMM-Newton , and Swift 
da ta is largely cons istent with the X-ray results reported 
by lAn ct ah ( 1201 3), but the current work provides im¬ 
provements and refinements. 

First, using longer Swift observations, we find that 
the orbital period of 1FGL J1018.6—5856 is 16.544 ± 
0.008 days, consistent with t he gamma-ra y measurement 
(16.531 ± 0.006 days; iColev et al. 1120141 ). When folded 
on the new period, the light curve shows two distinct 
features; a spike at phase 0 and a broad sinusoidal 
hump (Figure El), similar to those rep orted previously 
(lAckermann ct ahll2012t iAn et all 120131) . With the new 
period measurement, however, we find that the spike at 
phase 0 is a persistent feature and sh ows less orbit-to - 
orbit variability than was suggested bv lAn et al.l (12013 V 
Second, we clearly see the correlation betwe en flux and 
spectral hardness for which [An et al.l (120131) found only 
marginal evidence. This is possible thanks to more sensi¬ 
tive observations made with NuSTAR , Swift and XMM- 
Newton. 

Note that we combined all the Swift observations taken 
over a period of ~2000 days for the spectral analysis. If 
there is long-term (> P or b ) and/or short-term (10-100 ks) 
variability, the combined results may be incorrect. This 
is a concern because there are only a few observations 
per orbital phase bin, and individual exposure of the ob¬ 
servations is only ~ks. Furthermore, if the orbital period 
is not accurate or varies with time, phases of later ob¬ 
servations will change, introducing an additional error to 
the analysis. However, the agreement of the Swift mea¬ 
surements with the NuSTAR and XMM-Newton results 
suggests that the errors may not be large compared to 
the statistical uncertainties, having no significant impact 
on the results. For example, the significance of the corre¬ 
lation between flux and photon index is still >99% when 
adding 10% systematic uncertainty to the Swift measure¬ 
ments. 

The broadband X-ray spectra of 1FGL J1018.6—5856 


at phases 0 and 0.2 are well described with a 
power-law model in the 0 .5-40 keV band. Recently, 
IWaisberg fe Romani I (I2015T) find that, based on param¬ 
eter space consistent with radial velocity measurements, 
a neutron star model is preferred over a typical stellar 
mass black hole, although both classes are still allowed 
for 1FGL J1018.6—5856. We check to see if the source 
shows any evidence for accretion, such as line features 
or an exponential cutoff in the X-ray spectrum, as is 
often seen in neutron star HMXBs, and find none (Fig¬ 
ure a. Furthermore, we find no clear evidence for an 
exponential cutoff at -EVutoff < 70keV (for spectrum at 
phase 0). We, therefore, set the 90% lower limit for 
-^cutoff t o be 34-39 keV for e- folding energies of 6-12keV 
(Ef, see lCoburn et al. 1I2002L for the range of Ef of neu¬ 
tron star HMXBs). This lower limit is large for a neu¬ 
tron star HMX B (typical -E/utoff ~ 10 — 20keV; e.g., see 
iCoburn et al. l[2002h . Note that the X-ra y pulsar X Per 
also known as 4U 0352+309) for which ICoburn et al. I 
20021) did not find a cl ear spectral cutoff turn ed out to 
have a cutoff at 69 keV dLutovinov et al.ll2012f) , which is 
comparable to the energy under which we did not find 
any evidence for a spectral cutoff in 1FGL J1018.6—5856. 
Also, high cut off energies > 70 keV have been seen in black 
hole binaries dGrove et al.lIl999h . Therefore, we cannot 
clearly rule out the possibility that 1FGL J1018.6—5856 
is a black hole binary or a neutron star bianry with 
unusually high cutoff energy based only on the spec¬ 
tral cutoff. Nevertheless, the continuum spectrum of 
X Per or other X-r ay binaries is very complex (e.g., 
ICoburn et al. I l2002f) while we see a simple power-law 
spectrum for 1FGL J1018.6—5856. This suggests that 
1FGL J1018.6—5856 may be a non-accreting neutron 
star system, which has also bee n suggested f or another 
gamma-ray binary LS 5039 (e.g.. lTorresjl2(frni . 

In analogy to LS 5039, we may identify the loca¬ 
tion of the sinusoidal X-ray peak at </> ~ 0.4 (Fig¬ 
ure m as i nferior conjunction, and the gamma-ray peak 
at <t> ~ 0 (IColev et al. I2014f) as superio r conjunction 
dKishishita et al.l120091 : Abdo et, all I2009T) . Then, the 
phase difference of ~ 0.4 between the two conjunc¬ 
tions implies that the orbit is eccentric. We note, how¬ 
ever, that it is not clear whether the X-ray and gamma- 
ray peaks are physically related to the conjunctions or 
the apastron/periastron passages, and that alignment of 
the X-ray and gamma-ray peaks with inferior and supe¬ 
rior conjunctions may not be precise. Therefore, more 
observations and detailed modeling are required in order 
to draw a firm conclusion. 

We find that the X-ray spectral properties of 
1FGL J1018.6—5856 clearly show orbital modulation. 
Pulsar models for gamma-ray binaries often attribute 
such orbital modulation with orbital variation in the ad i- 
abatic cooling timescale dKhangulvan et afll2007l 2008 ). 
the electron injection spectr um, or the loca tion and the 
shape of the wind nebula dDubus I I2006T) . The pul¬ 
sar models have been applied to the similar system 
LS 50 39 (e.g., spectral variability and recu rring X-ray 
flares; iKishishita et al.l [20091 : IAn et, al.1120131 b and have 
reproduced the overall spectral energy distribution (e.g., 
iDubus 1 120061) . However, whether or not these mod¬ 
els can explain the spiky feature at phase 0 we see in 
1FGL J1018.6—5856 needs to be investigated. 

We note that the high-flux state observed with Swift 
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at phase 0 (Figure [5}i and b) is not reproduced in other 
observations taken at the same phase. It may be because 
the two observations in the high-flux state were made in 
a very narrow phase interval and the later observations 
did not cover that phase interval. In order to see if this is 
the case, we first verified that the high-flux state was not 
produced by short timescale variability (~ks); it lasted 
for the full duration of the exposures of the observations 
(24 ks at MJD 55585.7 and 7 ks at MJD 55618.7 for Obs. 
IDs 00031912004 and 00031912011, respectively) which 
cover a phase interval of A 4> = 0.022 (at (f> = 0.034 ± 
0.011 for P or b = 16.54days). We then measured the 
phases of the other observations. We find that there are 
four observations made at the high-flux phase interval, 
and none of them was in the high-flux state. Since the 
phase of an observation can change significantly for a 
different orbital period, we further varied P 0 rb within 
the measurement uncertainty of 0.01 day, and find the 
same result. This suggests that there is orbit-to-orbit 
flux variability at phase 0. 

We find that the duration of the high-flux state is 
longer than 24 ks and shorter than 1.8 days. The min¬ 
imum duration is set to be 24 ks because the high-flux 
states last during the observation (see above). The 
maximum duration is set to be the interval between a 
high-flux state and the next non-high-flux observation, 
wh ich is 1.8 days f or both high-flux states. As noted 
bv lAn et al.l (1201311 . the observational properties of the 
flare such as duration and o rbital repeatab ility look more 
like that of LS 5039 (e.g., IKishishita et all 2009f) than 
those of LS I +61°303 fe.g. JLi et al. 1 12011b!: I Smith et~aH 
l2009lf . This may support the idea that the flares 
are p roduced by dump i ness of the stel lar wind (e.g., 
iZdziarski. Ncronov fc Chernyakova 1120101) since the stel¬ 
lar companion (Be star) of LS I +61°303 is differ¬ 
ent from those (O stars) of 1FGL J1018.6—5856 and 
LS 5039 as the flare properties do. However, how 
the dumpiness produces flares at one orbital phase for 
1FGL J1018.6—5856 or LS 5039 but not at random or¬ 
bital phases needs to be further investigated. 


5. CONCLUSIONS 

We present results of NuSTAR , Swift, and XMM- 
Newton observations of the gamma-ray binary 
1FGL J1018.6—5856. Using the Swift data, we measured 
the orbital period of the source to be 16.544±0.008 days, 
in agreement w ith the refined gamma-ray measurement 
of IColev et al. I (12014T) . The new period is only slightly 
different from that used in our previous X-ray study, and 
hence our spectral and temporal analysis results agree 
well with the previous X-ray measurements. We find that 
the flux enhancement at phase 0 occurs more regularly 
in time than was suggested previously based on Swift 
data. The new NuSTAR and XMM-Newton data allow 
us to show clearly the correlation between X-ray flux 
and spectral hardness of 1FGL J1018.6—5856. Finally, 
the broadband X-ray spectrum of 1FGL J1018.6—5856 
suggests that it may not be an accretion-powered system. 
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